/* 
This is a very brief program in order to test the ewald code. 

In the header file 'ewald.h' there is a maximal limit for the 
reciprocal space cutoff, which (if required) has to be adjusted.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//
//
///* maximal precision (used in cutoff) */
//#define PREC 1.0e-30
//
///* precision of alpha */
//#define ALPHA_OPT_PREC 1.0e-10
//
//double rcut  = 4.9;      /* real space cutoff. Note that this must be smaller than L/2, if minimum image convention is to be used! */
//int    kmax  = 24;        /* maximal k-vector */
//double Ewaldalpha=1.0;    /*Value of the Ewald Splitting parameter */
//
///* number of particles: */
//#define NP 8
//
///* sum of square of charges */
//double Q2 = 8.0;
//
///* dielectric constant of 'surrounding' medium */
//double epsilon = 1.0;   /* = vacuum */
//
///* Define three arrays for the particle-coordinates: */
//double x[NP], y[NP], z[NP];
//
///* Define one array for the particle-charges: */
//double Q[NP];
///* Coulomb-energy of the box */
//double E_Coulomb;
//
///* Length of the system */
//double L = 10.0;
//
///* real space cutoff. Note that this must be smaller than L/2, if
//minimum image convention is to be used! */
//double rcut2;   /* rcut^2 */
//
//
///* prefactor of the Coulomb potential */
//double coulomb_prefactor = 1.0;
//
///* Berechnet das Quadrat von double x, gibt double zurueck. */
//double SQR(double x) { return x*x; }
//
///* Wert von double x runden, gibt double zurueck. */
//double dround(double x) { return floor(x+0.5); }
//
//#define     PI     3.14159265358979323846264338328
//#define   wupi     1.77245385090551602729816748334
//
///* constant to use when alpha is to be optimized */
#define ALPHA_OPT 0.0

#include "ewdClass.h"

int main(void)
{
	//int i;

	//double s1 = 0.50*L;
	//ewald *ewd=new ewald(6,8);

	///* using the formula of Benson-Mackenzie on gets for the Madelung
	//constant of the NaCl lattice: */
	//const double madelung_NaCl = -1.747564594633182190636212035544397403481;

	///* define the particle coordinates and charges such such that this
	//gives an NaCl-lattice upon periodic replication: */

	//x[ 0] = 0;  y[ 0] = 0;  z[ 0] = 0;  Q[ 0] =  1.0;
	//x[ 1] = s1; y[ 1] = s1; z[ 1] = 0;  Q[ 1] =  1.0;
	//x[ 2] = s1; y[ 2] = 0;  z[ 2] = s1; Q[ 2] =  1.0;
	//x[ 3] = 0;  y[ 3] = s1; z[ 3] = s1; Q[ 3] =  1.0;
	//x[ 4] = s1; y[ 4] = 0;  z[ 4] = 0;  Q[ 4] = -1.0;
	//x[ 5] = 0;  y[ 5] = s1; z[ 5] = 0;  Q[ 5] = -1.0;
	//x[ 6] = 0;  y[ 6] = 0;  z[ 6] = s1; Q[ 6] = -1.0;
	//x[ 7] = s1; y[ 7] = s1; z[ 7] = s1; Q[ 7] = -1.0;

	//rcut2 = SQR(rcut);

	///* initialize the Ewald-algorithm */
	//Ewaldalpha=ALPHA_OPT;


	//ewd.Ewald_init(L,NP,Ewaldalpha,rcut,kmax,coulomb_prefactor,Q2,epsilon);

	//printf("Tuned Ewald: alpha_opt=%e KPerr=%e\n", Ewald_get_alpha(), Ewald_KPerr());

	///* set forces and energy to 0 at the beginning: */
	//for (i = 0; i < NP; i++) Fx[i] = Fy[i] = Fz[i] = 0.0;
	//E_Coulomb = 0.0;   

	///* calculate the real space contribution to the coulomb energy/forces */
	//E_Coulomb += Ewald_r_space(x,y,z,Q,1,Fx,Fy,Fz);

	///* calculate the reciprocal space contribution to the coulomb energy/forces */
	//E_Coulomb += Ewald_k_space(x,y,z,Q,1,Fx,Fy,Fz);

	///*     WE DO NOT CALCULATE THE DIPOLE CONTRIBUTION !!!
	//...but if we wished, this would be done like:
	//E_Coulomb += Ewald_dipol(x,y,z,Q,1,Fx,Fy,Fz)
	//*/

	///* print the results. Note that in a perfect crystal the forces on
	//all particles should be identically zero! */

	//for (i = 0; i < NP; i++)
	//	/*printf("Fx[%d]=% le\tFy[%d]=% le\tFz[%d]=% le\n",i,Fx[i],i,Fy[i],i,Fz[i])*/;

	//printf("\nCoulomb-energy of the Box: %le\n",E_Coulomb);
	//printf("according to that, you predict a  Madelung constant = %14.10lf\n",L*E_Coulomb/NP);
	//printf("relative error, respect real Madelung constant: %le\n",(L*E_Coulomb/NP-madelung_NaCl)/madelung_NaCl);

	///*E_Coulomb += Ewald_cut();
	//printf("\n");
	//printf("relative error (with correction): %le\n",
	//(L*E_Coulomb/NP-madelung_NaCl)/madelung_NaCl);

	//*/
	//delete ewd;
	//return 0;
}







